A5. Estimate the model: Spatial smoothing with unknown changes in the pattern of response

Description

Fit a spatial semiparametric model based on splines including unknown changes in the pattern of response.

Usage

scp(formula, data, initial = NULL, contrasts = NULL,
        model = "exponential", fix.nugget = FALSE, fix.kappa = FALSE,
        nugget.tol = 1e-15, angle = 0, ratio = 1, use.reml = FALSE,
        use.profile = TRUE, chMaxiter = 20, control = list())

Arguments

formula

formula. An expression to specify the model to fit. See 2. Mean model.

data

sss class object. A dataset object generated by any of the commands as.sss or create.sss.

initial

named list. The starting values for the covariance parameters of the model. If initial=NULL then it is used an internal grid search to define the starting values.

contrasts

character. A contrast method for factor covariates. Default to contr.treatment.

model

character. Name of the semivariogram model to estimate for the spatial dependence. See Semivariogram Model.

fix.nugget

logical or numeric. If FALSE the nugget $\tau^2$ is estimated. If fix.nugget is a numeric value then the nugget $\tau^2$ is set to the value defined for fix.nugget.

fix.kappa

logical or numeric vector. If FALSE the parameters $\kappa$ are estimated. If fix.kappa is a numeric vector then $\kappa$ is set to the values of the vector defined for fix.kappa.

nugget.tol

numeric. Threshold for microscale spatial variations to define the nugget effect. Default to 1.0e-15. Do not modify unless know what is being doing.

angle

numeric. Angle for geometric anisotropy. Note that this overwrites any specification for aniso.angle in s2D.

ratio

numeric. Ratio between $[0,1]$ for geometric anisotropy. Note that this overwrites any specification for aniso.ratio in s2D.

use.reml

logical. For using REML estimation set to TRUE, for ML estimation set to FALSE (default).

use.profile

logical. For profiling set to TRUE (default).

chMaxiter

integer. Maximum number of iterations for the loop estimating changes in the pattern of response.

control

named list. Options to control the optimization. See argument control in command optim.

1. Semiparametric model

Assume that the response variable admit the trend surface model $$Y(s) = a^T b + g(s) + \epsilon(s)$$ where $a$ is a known vector of covariates and $b$ their coefficients; $g(s)$ is a deterministic bivariate spline and $\epsilon(s)$ is a Gaussian spatial process (GSP) with mean zero and covariance depending only on the distance $h$ and given by $Cov(\epsilon(s+h),\epsilon(s))$. This model is also called a trend surface model. Given $n$ observed locations $s_1,\ldots,s_n\in S\subset\Re^2$ in a two-dimensional space, then the model is $$Y = Ab + g + \epsilon$$ where $Y = (Y(s_1), \ldots, Y(s_n))^T$, $A$ is the known matrix of covariates, $g = (g(s_1),\ldots,g(s_n))^T$ and $\epsilon = (\epsilon(s_1),\ldots,\epsilon(s_n))^T$. The covariance matrix is given by $Cov(\epsilon,\epsilon) = \Sigma = \sigma^2 R$ with $R$ a valid correlation matrix. Thus $Y\sim N_n(\mu,\Sigma)$ where $\mu = Ab + g$ and the likelihood function is $L(b,g,\sigma^2,\theta;Y) = (2\pi)^{-n/2}|\Sigma|^{-1/2}\exp\{-\frac{1}{2}(Y-\mu)^T\Sigma^{-1}(Y-\mu)\}$ with $\theta$ the parameters that define the correlation matrix $R$.

2. Mean model

It can be defined by the commands:

linear: that defines the covariates in the matrix $A$. Note that more than one linear command can be defined. See linear.
cp: defines changes in the pattern of response by including the covariates $(z_d - \psi^{(0)}_d)\times 1\{z_d>\psi^{(0)}_d\}$ and $-1\{z_d>\psi^{(0)}_d\}$ for $d=1,\ldots,G$ into the matrix $A$. Note that more than one cp command can be defined. See cp.
s2D: define the bivariate spline $g$. Note that only one s2D command can be defined. See s2D.

3. Covariance model and nugget effect

Given a distance $h$ define $u = ||T^{1/2}_{\texttt{\tiny angle, ratio}}h|| = (h^T T_{\texttt{\tiny angle,ratio}} h)^{1/2}\in\Re$ where $T_{\texttt{\tiny angle, ratio}}$ is a rotation matrix for geometric anisotropy. The errors are given by the process $\epsilon(s) = \eta(s) + \xi(s)$ where $\xi$ is a GSP with mean zero and covariance $$\begin{array}{ll}Cov(\xi(s),\xi(s+h)) & = C_\xi(u;\phi,\kappa) \\ & = \sigma^2_0\rho_\xi(u;\phi,\kappa)\end{array}$$ with $\rho_\xi(u;\phi,\kappa)$ the correlation function; and $\eta$ is a nugget effect with covariance $$\begin{array}{ll}Cov(\eta(s),\eta(s+h)) & = C_\eta(u;\tau^2,\texttt{tol.nugget}) \\ & = \tau^2\rho_\eta(u;\texttt{tol.nugget})\end{array}$$ with correlation function $\rho_\eta(u;\texttt{tol.nugget}) = 1\{u<\texttt{tol.nugget}\}$. Therefore the covariance of the process $\epsilon$ is given by $$\begin{array}{ll} Cov(\epsilon(s),\epsilon(s+h)) & = C_\epsilon(u;\sigma^2,\theta,\texttt{tol.nugget}) \\ & = \sigma^2 \rho_\epsilon(u;\theta,\texttt{tol.nugget})\end{array}$$ with correlation function given by $$\rho_\epsilon(u;\theta,\texttt{tol.nugget}) = (1-\rho_*)\rho_\eta(u;\texttt{tol.nugget}) + \rho_*\rho_\xi(u;\phi,\kappa)$$ where $\theta = (\rho_*,\phi,\kappa)^T$ are the parameters with $\rho_* = \sigma^2_0/\sigma^2$, $\sigma^2 = \tau^2 + \sigma^2_0$, and tol.nugget is the argument that controls the largest distance at which micro-scale variations can affect the observed outcome. By default tol.nugget is set to 1.0e-15. The parameters $\phi,\kappa$ define the correlation function of the process $\xi$ with $\phi$ usually called the range parameter and $\kappa$ depending on the model selected. The semivariogram can be expressed as $$\gamma_\epsilon(u;\sigma^2,\theta,\texttt{tol.nugget}) = \sigma^2(1-\rho_\epsilon(u;\theta,\texttt{tol.nugget}))$$ where $\tau^2$ is the nugget effect, $\sigma^2$ is the sill, and $\sigma^2_0$ is the partial sill. Note that when angle = 0 and ratio = 1 the matrix $T_{\texttt{\tiny angle,ratio}}$ is an identity matrix and $u = h$ so the correlation $\rho_\epsilon(u;\theta,\texttt{tol.nugget})$ is isotropic. Use different values for angle and ratio to define a geometric anisotropic correlation function. Then the covariance matrix $\Sigma = \sigma^2 R$ where $R$ is the correlation matrix originated from $\rho_\epsilon(u;\theta,\texttt{tol.nugget})$. It is possible to define the argument model=name where name is one of the following: matern, powered.exponential, spherical, wave, exponential, gaussian, cubic, circular, gencauchy, cauchy, RMmatern, RMwhittle, RMgneiting, and RMnugget. For .semiVar one of matern, gaussian, exponential, power, cubic, penta.spherical, spherical, wave, sin.hole, pure.nugget and identity. By default the covariance model is set to exponential with angle=0 and ratio=1.

4. Penalized maximum likelihood estimation

Estimation can be performed by maximisation with respect to $b,g,\sigma^2,\theta$, and $\alpha$ of the penalized log likelihood $$\ell_p(b,g,\sigma^2,\theta,\alpha) = \log(L(b,g,\sigma^2,\theta;Y)) - \frac{1}{2\sigma^2} J_\alpha(g)$$ where $J_\alpha(g) = g^T Q_\alpha g$ is the penalty and $Q_\alpha$ is the roughness matrix.

5. Penalties

Depending on the type of spline assumed for $g$ the penalty is defined differently depending on the roughness matrix $Q_\alpha$ which is given by:

Tensor product spline.: Given $\tau_{1,1},\ldots,\tau_{1,K_1}$ and $\tau_{2,1},\ldots,\tau_{2,K_2}$ the design points in each coordinate then $$Q_\alpha = \alpha_1 I_{K_2}\otimes Q_1 + \alpha_2 Q_2\otimes I_{K_1}$$ where $Q_1$, $Q_2$ are unidimensional roughness matrices from the design points in each coordinate and $\alpha_1$, $\alpha_2$ are smoothing parameters in each direction.
Thin plate spline.: Given the $n$ locations, $Q_\alpha = \alpha E$ where $\alpha$ is the smoothing parameter and the $n\times n$ matrix $E$ has elements $E_{i,j} = \vartheta(||s_i - s_j||)$ for $i,j=1,\ldots,n$ where $$\vartheta(u) = \left\{\begin{array}{ll}\frac{1}{16\pi}\times u^2\log(u^2) & \textrm{ , $u>0$} \\ 0 & \textrm{ , otherwise.}\end{array}\right.$$

6. Mixed model representation

The spline can be written as $g = X\beta + Zr$ with $\beta$ and $r$ the coefficients and $X$ and $Z$ design matrices conveniently defined. Then for the observed responses the model can be expressed as a the mixed model $$Y = Ab + X\beta + Zr + \epsilon$$ where $r \sim Normal(0,I_V)$ with $V$ the number of columns in $Z$. Then, $Y\sim N_n(\mu_m, \Sigma)$ where $\mu_m = Ab + X\beta$ and $\Sigma = \sigma^2 R$; and $Y|r \sim N_n(\mu,V)$ where $\mu = Ab + X\beta + Zr$ and $V = ZZ^T + \Sigma$. Let us denote $\vartheta = (b,\beta,\sigma^2,\theta,\alpha)^T$, then the conditional log-likelihood of the model is given by $$\ell(\vartheta|r) \propto -\frac{1}{2}\left\{\log|\Sigma|+(Y-\mu)^T\Sigma^{-1}(Y-\mu)\right\}$$ and the marginal log-likelihood is given by $$\ell(\vartheta) \propto -\frac{1}{2}\left\{\log|V|+(Y-Ab-X\beta)^T V^{-1}(Y-Ab-X\beta)\right\}.$$